In [1]:
data_dir = r'C:\Data\Antonio\Philip\081114Patch clamp\nanorods 630\fov2 - good\122055_take2 100Hz\\'
In [2]:
from __future__ import division
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
sns.set_context(rc={'lines.markeredgewidth': 1}) # workaround for a bug in matplotlib 1.4.2
In [3]:
from patchclamp import PatchDataset
import timetraces as tt
In [4]:
data = PatchDataset(data_dir)
The attributes data.voltage/current/time are resampled at the camera frame rate:
In [5]:
plt.plot(data.time[:100], data.voltage[:100], label='Voltage')
plt.ylabel('Voltage (V)')
plt.twinx()
plt.plot(data.time[:100], data.current[:100], label='Current',
color=sns.color_palette()[1])
plt.ylabel('Current (pA)')
plt.grid(False)
plt.xlabel('Time (s)')
Out[5]:
In [6]:
sns.set_style('dark')
In [7]:
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(data.video.mean(0), cmap='cubehelix')
plt.colorbar(im);
ax.set_title('Mean frame')
Out[7]:
In [8]:
import scipy.ndimage as ndi
In [9]:
mvideo = data.video.mean(0)
In [10]:
data.video.shape, data.video.dtype
Out[10]:
In [11]:
smooth_mvideo = ndi.gaussian_filter(mvideo, sigma=20)
smooth_mvideo.shape, smooth_mvideo.dtype
Out[11]:
In [12]:
smooth_frame = ndi.gaussian_filter(data.video[0].astype('float64'), sigma=20)
smooth_frame.shape, smooth_frame.dtype
Out[12]:
In [13]:
smooth_video = ndi.gaussian_filter(data.video.astype('float64'), sigma=20)
smooth_video.shape, smooth_video.dtype
Out[13]:
In [14]:
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(smooth_mvideo, cmap='cubehelix', vmin=104, vmax=109)
plt.colorbar(im);
ax.set_title('Mean background frame');
In [15]:
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(smooth_frame, cmap='cubehelix', vmin=104, vmax=109)
plt.colorbar(im);
ax.set_title('Background frame 0');
In [16]:
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(smooth_video.mean(0), cmap='cubehelix', vmin=104, vmax=109)
plt.colorbar(im);
In [17]:
corr_video = data.video.astype(np.float64) - smooth_video
corr_video.shape, corr_video.dtype
Out[17]:
In [18]:
nframe = 5
fig, ax = plt.subplots(figsize=(10, 6))
im = ax.imshow(corr_video[nframe], cmap='cubehelix', vmin=-50, vmax=60)
plt.colorbar(im);
ax.set_title('Background-corrected frame %d' % nframe);
In [20]:
%matplotlib qt
In [21]:
from figscroller import FrameScroller
Explore the video frame by frame:
In [22]:
fig, ax = plt.subplots(figsize=(14, 9))
im = ax.imshow(corr_video[0], cmap='cubehelix')
plt.colorbar(im)
scroller = FrameScroller(fig, im, corr_video)
points = plt.ginput(10, timeout=0)
points
Out[22]:
Or use only the mean frame:
In [12]:
#frame = corr_video.mean(0)
frame = data.video.mean(0)
fig, ax = plt.subplots(figsize=(14, 9))
im = ax.imshow(frame, cmap='cubehelix')
plt.colorbar(im)
points = plt.ginput(10, timeout=0)
points
Out[12]:
QD on patched cell:
In [8]:
points_patched = \
[(92.954301075268802, 122.91935483870967),
(93.411290322580612, 127.26075268817203),
(91.81182795698922, 132.05913978494624),
(128.14247311827953, 136.17204301075267),
(129.28494623655911, 135.02956989247309),
(139.1102150537634, 125.2043010752688),
(147.10752688172039, 66.938172043010752),
(101.4086021505376, 88.188172043010752),
(100.72311827956986, 89.787634408602145),
(154.64784946236554, 95.956989247311824)]
points_patched = \
[(93.182795698924707, 123.14784946236558),
(93.639784946236517, 127.03225806451611),
(92.040322580645125, 132.05913978494624),
(128.82795698924727, 135.71505376344084),
(147.33602150537629, 66.938172043010752),
(101.4086021505376, 88.188172043010752),
(100.95161290322577, 89.55913978494624),
(139.1102150537634, 124.97580645161288),]
QD unpatched:
In [9]:
points_unpatched = \
[(161.04569892473114, 109.89516129032256),
(163.78763440860212, 112.40860215053762),
(226.39516129032253, 41.118279569892479),
(246.95967741935479, 133.6586021505376),
(26.233870967741929, 140.0564516129032),
(103.92204301075265, 48.887096774193537),
(16.86559139784945, 58.712365591397855),
(102.32258064516125, 3.6451612903225907),
(104.15053763440858, 163.13440860215053),
(185.95161290322577, 148.96774193548384)]
points_unpatched = \
[(226.68894009216589, 41.434907834101381),
(193.39400921658986, 139.75288018433179),
(185.95161290322582, 149.54550691244236),
(247.05760368663599, 133.48559907834101),
(26.135944700460826, 140.53629032258061),
(6.5506912442396299, 178.53168202764974),
(104.08525345622121, 48.877304147465452),
(16.735023041474648, 58.278225806451616),
(38.670506912442391, 59.06163594470047),
(218.07142857142858, 5.0063364055299644)]
In [10]:
%matplotlib inline
In [11]:
sns.set_style('dark')
In [12]:
fig, ax = plt.subplots(figsize=(10, 6))
for point in points_patched:
plt.plot(point[0], point[1], 'r+')
im = ax.imshow(data.video.mean(0), cmap='cubehelix')
plt.colorbar(im);
ax.set_title('QD on patched cell')
fig, ax = plt.subplots(figsize=(10, 6))
for point in points_unpatched:
plt.plot(point[0], point[1], 'r+')
im = ax.imshow(data.video.mean(0), cmap='cubehelix')
plt.colorbar(im);
ax.set_title('Unpatched QD');
In [24]:
import timetraces as tt
In [79]:
%matplotlib inline
In [13]:
from matplotlib import mlab
In [16]:
nfft = 256
noverlap = 128
cycle_mean = False
clip_radius = 2
In [17]:
timetrace = tt.get_timetrace_circle(data.video, points_patched[1],
clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
In [18]:
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
chunked_signal.shape
Out[18]:
In [19]:
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
fft.shape
Out[19]:
In [20]:
frequency = np.arange(fft.shape[0])*max_freq/(fft.shape[0]-1)
frequency.shape
Out[20]:
In [21]:
# Pixel aspect ratio for the spectrogram
aspect = 0.13*(fft.shape[0]/(fft.shape[1]))
aspect
Out[21]:
If there is signal it will show up in the real part of the FFT, since the signal is always in phase and starts with an ON (or OFF) period.
In [22]:
spec = np.abs(fft.real)**2
fdrop = 4
fig, ax = plt.subplots(2, 1, figsize=(18, 4))
freqbin = frequency[1] - frequency[0]
ax[0].imshow(spec[fdrop:].T, interpolation='none', cmap='cubehelix', aspect=aspect,
extent=(-freqbin/2 + frequency[fdrop], frequency[-1] + freqbin/2, -0.5, spec.shape[1] + 0.5))
ax[0].tick_params(length=5, direction='out')
ax[0].grid(False)
ax[1].plot(frequency[fdrop:], spec[fdrop:].mean(1), marker='.')
ax[1].set_xlabel('Frequency (Hz)');
ax[1].set_xlim(frequency[fdrop], frequency[-1])
Out[22]:
In [23]:
video = data.video
cycle_mean = False
clip_radius = 2 # For circular selection
nfft = 256
noverlap = 128
fdrop = 4
fig, ax = plt.subplots(len(points_patched), 1, figsize=(18, 2*len(points_patched)),
sharex=False, sharey=False)
fig.subplots_adjust(wspace=0.04, hspace=0.03)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/8
im = ax[i].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1] + 0.5),
aspect=aspect)
ax[i].tick_params(direction='out', length=5)
ax[i].grid(False)
ax[-1].set_xlabel('Frequency (Hz)')
ax[0].set_title('Spectrogram');
In [24]:
video = data.video
cycle_mean = True
clip_radius = 2 # For circular selection
nfft = 256
noverlap = 128
fdrop = 32
fig, ax = plt.subplots(len(points_patched), 1, figsize=(18, 2*len(points_patched)),
sharex=False, sharey=False)
fig.subplots_adjust(wspace=0.04, hspace=0.03)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/16
im = ax[i].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1] + 0.5),
aspect=aspect)
ax[i].tick_params(direction='out', length=5)
ax[i].grid(False)
ax[-1].set_xlabel('Frequency (Hz)')
ax[0].set_title('Spectrogram');
In [25]:
video = data.video
cycle_mean = True
clip_radius = 2 # For circular selection
nfft = 32
noverlap = 24
fdrop = 4
fig, ax = plt.subplots(len(points_patched), 1, figsize=(18, 2*len(points_patched)),
sharex=False, sharey=False)
fig.subplots_adjust(wspace=0.04, hspace=0.03)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/2
im = ax[i].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1] + 0.5),
aspect=aspect)
ax[i].tick_params(direction='out', length=5)
ax[i].grid(False)
ax[-1].set_xlabel('Frequency (Hz)')
ax[0].set_title('Spectrogram');
In [27]:
sns.set_style('darkgrid')
In [29]:
video = data.video
cycle_mean = True
clip_radius = 1.5 # For circular selection
nfft = 16
noverlap = 8
fdrop = 1
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/0.4
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [30]:
video = data.video
cycle_mean = False
clip_radius = 2 # For circular selection
nfft = 16
noverlap = 8
fdrop = 1
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/0.4
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [31]:
video = data.video
cycle_mean = False
clip_radius = 2 # For circular selection
nfft = 32
noverlap = 8
fdrop = 1
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/0.4
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [32]:
video = data.video
cycle_mean = False
clip_radius = 2 # For circular selection
nfft = 64
noverlap = 8
fdrop = 4
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/1
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [33]:
f0 = 100
sim_signal = 1 + np.cos(2*np.pi*f0*data.time - np.pi/4)
In [34]:
plt.plot(data.time[:10], sim_signal[:10], '-o')
Out[34]:
In [37]:
video = data.video
cycle_mean = True
clip_radius = 2 # For circular selection
nfft = 64
noverlap = nfft - 8
fdrop = 9
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
# Add simulated signal
timetrace += 1*sim_signal
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/4
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
#ax[i, 1].set_ylim(0, 1e6);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [40]:
video = data.video
cycle_mean = True
clip_radius = 2 # For circular selection
nfft = 256
noverlap = 128
fdrop = 32
fig, ax = plt.subplots(len(points_patched), 2, figsize=(20, 2*len(points_patched)))
fig.subplots_adjust(wspace=0.06)
for i, point in enumerate(points_patched):
timetrace = tt.get_timetrace_circle(video, point, clip_radius=clip_radius)
# Add simulated signal
timetrace += 1*sim_signal
max_freq = data.camera_rate*0.5 # Nyquist frequency
if cycle_mean:
timetrace = timetrace.reshape(timetrace.size/2, 2).mean(1)
max_freq /= 2
chunked_signal = mlab.stride_windows(timetrace, nfft, noverlap=noverlap)
fft = np.fft.rfft(chunked_signal, axis=0)
fft[1:-1] *= 2
spectr = np.abs(fft.real)**2
freq = np.arange(spectr.shape[0])*max_freq/(spectr.shape[0]-1)
dfreq = freq[1] - freq[0]
aspect = spectr.shape[0]/spectr.shape[1]/8
im = ax[i, 0].imshow(spectr[fdrop:].T, cmap='cubehelix', interpolation='none',
extent=(freq[fdrop] - 0.5*dfreq, freq[-1] + 0.5*dfreq, -0.5, spectr.shape[1]+0.5),
aspect=aspect)
ax[i, 0].tick_params(direction='out', length=5, right=False, left=False)
ax[i, 0].grid(False)
ax[i, 1].plot(freq[fdrop:], spectr[fdrop:].mean(1), marker='s')
ax[i, 1].axvline(100, ls='--', color=sns.color_palette()[1])
#ax[i, 0].set_xticks(freq[::4]);
ax[-1, 0].set_xlabel('Frequency (Hz)')
ax[-1, 1].set_xlabel('Frequency (Hz)');
ax[0, 0].set_title('Spectrogram')
ax[0, 1].set_title('Mean spectrum');
In [ ]: